voles

Preperations

library(volesModel)

## This is "synlik"  0.1.2

library(plotly)
library(matrixStats)


data(voles_data)
voles_data$date<-voles_data$year+0.25
voles_data[which(voles_data$season=="autumn"),"date"]<-voles_data[which(voles_data$season=="autumn"),"date"]+0.5
voles_data<-voles_data[order(voles_data$date),]

Plot our data

ax<- list(
  title = "",
  # titlefont = f1,
  showticklabels = TRUE,
  tickangle = 0
  # tickfont = f2,
  # exponentformat = "E"
)

ay<- list(
  title = "Voles Trapping Index",
  titlefont = "Open Sans",
  showticklabels = TRUE,
  tickangle = 270,
  # range = c(0.1, 0.6),
  tickfont = "Times New Roman"
  # exponentformat = "E"
)

p <- plot_ly(data = voles_data, x = ~date, y = ~popDensity,type = 'scatter', mode = 'lines+markers') %>%
  layout(xaxis = ax, yaxis = ay, showlegend = FALSE)
p

Functions

Summary Statistics

volesStats  <- function(x, obsData, ...){

  

  if (!is.matrix(x)) x <- matrix(x, 1, length(x))
  tx <- t(x)

  X0 <- t(orderDist(tx, obsData, np=3, diff=1)) ## cubic regs coeff of ordered diffs
  X0 <- cbind(X0, t(nlar(t(x),lag=c(6, 6, 6, 1, 1), power=c(1, 2, 3, 1, 2)))) ## least square coeff
  X0 <- cbind(X0, rowMeans(x)) # mean values of Y, # of 0's
  X0 <- cbind(X0, t(slAcf(tx, max.lag = 5)))  ## autocovariances up to lag 5 (the first element is the variance)
  X0 <- cbind(X0, apply( t( abs( apply( sign( apply(x,1,diff) ), 2, diff ) )), 1, sum)/2) #Number of turning points
  X0 <- cbind(X0, rowMedians(x)) # Mean-median
  
  
  
  r=character()
  for (i in 1:3) r[i]<-paste("cubic_",i,sep = "")
  for (i in 4:8) r[i]<-paste("ols_",i-3,sep = "")
  r[9]<-"mean"
  r[10]<-"variance"
  for (i in 11:15) r[i]<-paste("autocov_",i-10,sep = "")
  r[16]<-"turning"
  r[17]<-"mean_med"
  
  colnames(X0)<-r
  return(X0)
} 

Simulator

Priors:

voleFullPrior_sl <- function(input){ sum( input ) +
                               dnorm (exp(input[1]), 5, 1, log = TRUE)          +
                               dnorm (exp(input[2]), 1, 1, log = TRUE)          +
                               dexp  (exp(input[3]), 7, log = TRUE)             +
                               dgamma(exp(input[4]), 4, 40, log = TRUE)         +
                               dnorm (exp(input[5]), 15, 15, log = TRUE)        +
                               dnorm (exp(input[6]), 0.04, 0.04, log = TRUE)    +
                               dnorm (exp(input[7]), 1.25, 0.5, log = TRUE)     +
                               dunif (exp(input[8]), 0.5, 30, log = TRUE)       +
                               dunif (exp(input[9]), 20, 300, log = TRUE)       }

How does simulation work?
1.) Set Starting values for weasle and vole abundance (e.g. 0)
2.) Sample from the prior the parameters.
3.) Create a whole dataset (T=90) with this priors and the starting values;
update weasle and vole abundance at each t
4.) Calculate summary statistics and accept or reject parameters according to real data summary statistics.
5.) Do steps 2.) to 4.) M times.
6.) For plotting: Use posterior mean of \(\theta\) to draw from \(Pois(\theta n_t)\) to create \(Y_t\)